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Techniques  useful  for  efficiently  time-stepping  Galerkin  methods  for  various 

types  of  time-dependent  partial  differential  equations  are  presented  and  analyzed. 

Second-order  quasilinear  hyperbolic  problems  with  smooth  solutions  are  studied  as 

a  simple  model  problem  for  illustrating  the  widely  applicable  techniques.  The 

procedure  involves  the  use  of  a  preconditioned  iterative  method  for  approximately 

solving  the  different  linear  systems  of  equations  arising  at  each  time-step  in  a 
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discrete-time  Galerkin  method.  Optimal  order  L  spatial  errors  and  almost  optimal 
order  work  estimates  are  obtained  for  the  second-order  hyperbolic  equation. 
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SIGNIFICANCE  AND  EXPLANATION 


Recently  very  efficient  procedures  have  been  developed  for  time-stepping  non¬ 
linear  partial  differential  equations  of  parabolic  type  by  the  author  and  others. 

This  report  extends  some  of  these  techniques  to  various  types  of  partial  differential 
equations  which  are  second-order  in  the  time  derivative.  Equations  to  which  these 
techniques  can  be  applied  have  been  used  as  models  for  vibrational  problems#  for 
dynamics  of  rotating  fluids,  for  nonlinear  viscoelasticity,  and  for  other  physical 
problems . 

Significant  amounts  of  computation  are  saved  by  an  iterative  time-stepping 

procedure  for  approximating  the  solution  of  the  large  systems  of  linear  equations 

produced  by  a  Galerkin-type  numerical  procedure.  Instead  of  factoring  a  different 

large  matrix  at  each  time  step  to  solve  the  linear  equations  exactly,  a  new  matrix 

is  factored  periodically  and  used  as  a  preconditioner  in  an  iterative  procedure. 

Work  estimates  that  are  presented  indicate  how  often  the  preconditioner  must  be 

updated  to  obtain  near-optimal  amounts  of  work.  Very  few  iterations  are  then 

required  at  each  time  step  since  the  iterative  procedure  is  just  a  stabilizing 

process  for  the  underlying  time-stepping  procedure.  A  complete  error  analysis  is 

presented.  _  _ 
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ON  EFFICIENT  TIME-STEPPING  METHODS  FOR  NONLINEAR 
SECOND  ORDER  HYPERBOLIC  PARTIAL  DIFFERENTIAL  EQUATIONS 


Richard  E.  Ewing^ 


1.  Introduction 

Ve  shall  consider,  as  a  model  problem,  efficient  procedures  for  time- 
stepping  Galerkin  methods  for  approximating  smooth  solutions  of  quasilinear 
second-order  hyperbolic  equations.  The  techniques  presented  in  this  analysis 
can  be  applied  to  various  generalized  wave  equations  which  are  used  as  model 
equations  for  many  different  types  of  vibrational  problems.  We  consider  the 
problem  of  approximating  the  smooth  solution  u  “  u(x,t)  which  satisfies 


a)  c(x)  — r-  V-Ia(x,u)VuJ  *  f(x,t,u), 
3t 


X  e  0,  t  €  J  , 


(1.1) 


b)  u(x,0)  -  Uq(x)  , 


c)  (x,0)  »  Vq(x)  , 


X  e  Q  , 


X  c  0  « 


d)  a(x,u)  ~  -  g(x,t)  , 


X  e  an,  t  e  J  , 


where  0  is  a  bounded  domain  in  IR  ,  d  3,  with  boundary  Jfl,  v  is  the 
outward  unit  normal  to  9ft,  J  =  (0,TJ,  and  c,  a,  f,  u^,  v^,  and  g  are 
prescribed.  Ke  shall  first  present  a  Crank-Nicolson-Galerkin  approximation 
to  (1.1)  which  produces  a  different  linear  system  of  equations  to  be  solved 
at  each  time  step.  Procedures  of  this  type  have  been  analyzed  in  [7,8,11,3]. 
Our  modification  of  the  basic  procedure  will  consist  of  using  a  preconditioned 
iterative  procedure  for  only  approximating  the  solution  of  these  linear 
equations  at  each  ti'  j  step.  The  use  of  a  preconditioning  matrix  eliminates 
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the  need  to  refactor  a  new  matrix  at  each  time  seep,  while  the  iterative 
procedure  is  used  to  stabilize  the  resulting  algorithm.  Using  this  modifica¬ 
tion,  we  obtain  the  same  order  error  estimates  as  for  the  base  scheme  with 
greatly  reduced  computational  complexity.  We  obtain  very  nearly  optimal 
possible  work  estimates  for  our  procedure. 

The  techniques  presented  here  can  also  be  used  to  analyze  approxima¬ 
tion  procedures  for  initial-boundary  value  problems  for  equations  of  the  form 
-2  _ 

(1.2)  c(x)  — V*  (a(x,Vu)7u  +  b(x,Vu)  V 1  «=  f(x,t,u,Vu),  x  «  0»  t«  J  , 

with  appropriate  initial  and  boundary  conditions.  Equations  of  this  type 
have  been  used  as  models  in  nonlinear  viscoelasticity  and  hydrodynamics. 
Existence,  uniqueness,  and  stability  of  equations  of  this  type  have  been 
studied  by  Dafermos,  Greenberg,  MacCamy,  Hizel,  Showalter  and  others  (see 
[5,6,18,22,25]}.  The  coefficient  a  can  be  allowed  to  degenerate  to  zero 
In  (1.2). 

He  can  also  treat  approximations  of  solutions  of  equations  of  the  form 


(1.3) 


2  2 
c(x)  ^  -  V*fa(x,7u)Vu  +  b(x,Vu)7  —  +  e(x,7u)7  -  f(x,t,u,7u) 

at  I*  af'J 


X  c  Q,  t  «  J 


with  appropriate  initial  and  boundary  conditions.  Equations  of  this  type 
have  been  used  as  classical  vibration  models  -  (21 ,  §278]  *and  in  the  dynamics 
rotating  fluids  (20,23].  The  coefficients  a  or  b  can  be  allowed 
to  degenerate  to  zero  in  (1.3). 

Efficient  time-stepping  procedures  of  the  type  presented  here  have 
been  used  by  the  author  and  others,  for  pseudoparabolic  equations  in  (13], 
for  parabolic  equations  in  [10,4]  ,  and  for  systems  of  equations  used  to 
model  miscible  displacement  in  porous  media  in  [14,15,16]. 


In  Section  2  we  introduce  finite  element  spaces,  present  the  hypotheses 
on  (1.1)  and  its  solution  u,  discuss  an  elliptic  projection  of  u,  and 
present  various  Crank-Nicolson-Galerkin  methods  for  {1.1)-(1.3).  In 
Section  3  we  present  our  preconditioned  modification  of  the  base  method 
and  discuss  the  effect  of  the  iterative  stabilization  on  a  single  time  step. 
We  obtain  global  error  estimates  for  both  the  base  scheme  and  the  iterative 
modification  in  Section  4.  Section  5  contains  a  brief  description  of 
estimates  of  the  computational  complexity  of  the  methods  presented  in  the 
paper. 


2.  Preliminaries  and  Description  of  Galerkin  Methods 

l«*t  “  /  ^^dx,  =  /  ^\|ids,  and 

0  30 

I«et  W_(0)  be  the  Sobolev  space  on  0  with  norm 

8 


f,  ? 

l  loli' 


<k  "  3x 


L®(n) 


vlth  the  usual  modification  for  s  ■  When  s  ■*  2,  let 
ll'J'll  k  “  Ikll  ^  “  Iklljj-  If  VF  «  (F^^F^),  write  1Ipf||  ^  in  place  of 

\q  H  w^ 

1 

k  k'  *  H^{30)  denote  the  corresponding  Sobolev  space 

w  w 


on  30  with  norm  ||i;;| 


hN3»,  '  '*'»  ^  '♦'o’- 


Let  be  a  family  of  finite-dimensional  subspaces  of  (0)  with 

the  following  property: 


For  p  =  2  or  p  =  ”,  there  exist  an  integer  r  >  2 
and  a  constant  such  that,  for  1  ^  q  £  f  and 

(2.1)  inf  xll  0  +  Xll  ,}  IK h*’  . 

n  p  P  P 

We  also  assume  that  the  family  {A(^}  satisfies  the  following  so-called 
"inverse  hypotheses";  if 

_d 

*)  lull  .  1  KqH  ^  M  , 

(2.2)  ^ 

w  lUlli  1  • 

Q  as  follows  (with  (S)  denoting  the  collection  of  restrictions) 

1)  The  Neumann  problem  for  -A  +  l  on  Q  is  H  -regular. 

(S) 

2)  3Q  is  Lipschitz. 

Assume  the  following  regularity  for  c,  a,  a,  b,  e,  f  and  u: 

(Q)  1.  niere  exist  uniform  constants  such  that 

s)  0  <  a^  £  a(x,u)  i  a*  , 

b)  0  £  a*  <,  a(x,Vu)  i  a*  £  Kj^  , 

c)  0  <  £  c{x)  £  c*  £  , 

(2.3)  d)  0  <  <  b(x,Vu)  <  , 

•)  ®  £  b^  £  b(x,7u)  _<  , 

£)  0  <  e^  £  e(x,Vu)  £  Kj^  , 

g)  |f{x,t,u) I  £  . 

2.  The  functions  a  =  a(x,u),  a  *=  a(x,qj^,q2)*  b  =  b(x,qj,q,), 
b  =  b(x,qj^,q2),  e  =  e(x,qj^,q2)  ,  and  f  =  f{x,t,u)  are 
continuously  differentiable  with  respect  to  u  (respectively 
Vu)  and  have  a  uniform  bound,  K^,  satisfying  (for  i  ”  1,2) 


-4- 


(2.4) 


|al,  |al,  lb|,  |bl.  Id,  hi,  |fl,  Ifjl.  I^I.  1^1, 


.2- 

3  a  I 


1^1.  igi.  ifil. 

3u 


3qi 


Define 


(2.5) 


W’((a,b);X) 

P 


=  II  lh(-.t)|Ll 


X  **  Q 

*  W’{a,b) 


1  1  p»q  1 


Let  u,  the  solution  of  (1.1)  satisfy  the  following  regularity  assurcptions : 

3ul) 


R: 


1  2  r 
l'  (JjH  ) 


a^u 


3t 


L^(J;H*^) 


(2.6) 


b)  Hull .  3  + 

L  (J;H^) 


3u 

3t 


3^u 


L  (J;H  ) 


3t 


<  Kj  , 


L  (J;H  ) 


c) 


3^ 


3t^ 


1 

.3 

II 

^  1 

3  u 

▼  1 

3 

*  1  4 

3t 

,  ,  "at 

1^2  • 


2  2 
L  (J.-L'^) 


L  (J;L  )  L  (J;H  ) 

Similar  regularity  assumptions  must  be  satisfied  by  the  solutions  of  (1.2) 
and  (1.3)  but  we  shall  not  make  these  explicit  here. 

As  in  [261 ,  we  shall  introduce  an  auxiliary  elliptic  problem  to  aid  in 
our  analysis.  Define  W  in  to  be  the  unique  function  which,  for 


t  f  [0,T] ,  satisfies 

(2.7)  (a(u)VW,Vx)  +  (W,x)  «  (a(u)7u,7x)  +  (u,x) ,  X  < 

Shen  as  in  [26,7,12]  we  obtain  the  following  lemma. 


Lemma  2.1.  There  exists  a  constant  X3  «  K3 (n,a^ ,Kq,K^ .K^)  such 
that  for  2  <,  q  i  r,  n  *=  u  -  W,  and  s  «=  0  or  1, 


n 


-5- 


0  llnll  .  ^  1  Kjh’-'lloll  . 

L  (J;H  )  L 


b) 


in 

at 


(2.8) 


l!n 

3t^ 


1  ^3^ 


q-s 


{llull  3  , 

L  (J;h' 


) 


5u 

3t 


a^u 

at^ 


}■ 


We  also  make  the  assumption  on  **  that  there  exists  a 


constant  K.  such  that 

4 


(2-«l  I|W|I  .  .  *  I|VW||  .  .  .  |||2|| 


L  (JjL  ) 


L  (J;!.  ) 


at' 


*  ll'’l?ll 


IK,  • 


L  (J;L  )  L  (J;L  ) 

Sufficient  conditions  for  (2.9)  to  hold  can  be  found  in  [10,26)^  Also  as  in 
17,10,12]  we  can  obtain  the  following  lemma. 

Lemma  2.2.  There  exists  a  constant  “  K^({l,a^tKQ,K^,K^)  such  that 


(2.10) 


II  5^1 

+  II  — 

ii  2 

^  I  3 

» at  ' 

,  “  at 

1  Kg  • 


L  L  (J;H  ) 

We  shall  consider  discrete~time  Galerkin  approximations.  Let  At  >  0, 

M  ■  T/At  e  22,  and  t°  *  oAt,  o  «  ».  Also  let  =  (»"(x)  5  i)»(x,t"),  and 

n 


(2.11) 


..  .n+1  -,n  ^  ,n-l 

w  ■ 

*  »«■* 


We  shall  consider  Crank-Nicolson-Ualerkin  methods  for  our  base  time- 


stepping  procedure  for  each  of  the  equations.  Let  U  :  tbQ,...,t^} 


be  an  approx ir.at ion  to  the  solution  of  (1.1).  Assurdng  that  U  are 
known  for  k  ^  n,  determine  ^  by 

(2.12)  (cd^u'',x)  +  |a{u")V  - y! - ,7xj  =  (f  (t",u")  ,x)  +  <9  (t")  ,x>,  X'  \ 

Similarly,  we  define  our  approximation  to  the  solution  of  (1.2)  by 


(cd^U  ,x)  +  a(Vu")V  - - - -  .Vx  +  b(yu")V  - - - ,  Vx 

\  t  j  J  \  * 


,n+l  ..n-1 


(2.13) 


=  (f  (t",u'',Vu"),x)  +  <5(t''),x>,  Xf 


and  our  approximation  to  the  solution  of  (1.3)  by 


_  /  f„n+l  ^  r  j.  ^ 

(cd^u",x)  +  |a(7u'')v|^^ - ^ - ^],7x]  +  |b(yu")V  - — - ,  Vx] 

+  (e(Vu")VdJu",Vx)  =  (f(t",u".vu"),x)  +  <g(t").x>» 


3.  Iterative  Procedures 

In  this  section,  we  shall  present  the  linear  equations  arising  from 
(2.12)- (2.14) .  We  note  that  in  each  case,  the  coefficient  matrices  change 
with  each  time  step.  In  order  to  avoid  factorization  of  different  matrices 
at  each  time  step  for  the  solution  of  the  linear  equations,  we  shall  discuss 
an  iterative  method  for  approximating  their  solution.  The  analysis  pre¬ 
sented  here  will  extend  the  analysis  of  [10,13]  to  the  equations  (1.1)- (1.3) . 

Let  {y.}**  be  a  basis  for  M.  .  Let  0°  from  (2.12)  be  written  as 

i  i-1 


(3.1) 


^ -  I  • 

i-1  * 


Ve  then  see  that  using  (3.1) ,  (2.12)  can  be  %nritten  as 


"(C)  (C 


-  c")  -  C(c"  -  c"‘^) 


(3.2) 


-  a"{C)(c"  +  c"~^)  +  (At)^F"{C)  =  Rj{C) 


7 


4 


where 


a)  C  •=  ((cp^,p^))  » 


M 


(3.3)  b)  A"(t)  •=  ((a(  I  c"pJVy  ..Vu.)) , 

t=l  ^ 


and 


M 


c)  f"(C)  =  ((f(t",  I  C?v,).Pi)  +  <g(t"),y.>)  . 
1  1=1 


Similarly,  (2.13)  can  be  written  as 


(3.4) 


[■ 


C  +  a"(C)  + 


~  b"(C)]  (C 

-  [c-f  - 


-  C'') 


a”(5)  (t"  +  c""^)  + 


and  (2.14)  can  be  vrritten  as 


(3.5) 


[c  ♦  £"(5)  +  ■^^a"{C)  +  ~  B"(5)j  (C 

-  +  e"(5)  -  —  B"(C)j  (5”-  C"’^)  - 


n+1  _ 


a" (5)  (t"  +  5"“^)  + 


where  b"  and  e"  are  defined  as  in  (3.3.b)  with  the  coefficient  a 
replaced  by  b  and  e  respectively  and  F^  i*  defined  in  an  analogous 
Bianner  to  F^. 

Mote  that  since  the  matrices  A*,  B*,  and  E**  change  with  time, 
straightforward  solution  of  (3.2),  (3.4),  or  (3.5)  would  involve  the  factoriza¬ 
tion  of  new  matrices  at  each  time  step.  Instead  of  solving  (3.2)  exactly, 
tre  shall  approximate  the  solution  by  using  an  iterative  procedure  which 
has  been  preconditioned  by 


(3.6) 


»®  _  r  4.  tPtry 

L  ■  C  — - —  A  (5)  . 


Similarly,  for  (3.4)  and  (3.5)  we  shall  precondition  with 

£®  -  C  ♦  A®(C)  +  Y 


-8- 


L°  =  C  +  E°(^)  +  +  Y  , 

respectively.  The  preconditioning  process  eliminates  the  need  for  factoring 
new  matrices  at  each  time  step,  while  the  iterative  procedure  stabilizes  the 
resulting  proble-m.  The  st^d)ili2ation  process  requires  iteration  only 
until  a  predetermined  norm  reduction  is  achieved. 

Let  the  approximation  of  u”  from  (2.12)  produced  by  only  appro .ximateiy 
solving  (3.2)  using  the  preconditioner  (3.6)  be  denoted  by 


(3.7) 


v"-  I 

i-1 


A  Starting  procedure  for  determining  and  will  be  discussed  later. 

Assuming  that  these  quantities  are  known,  we  shall  determine  n  ^  1, 

using  a  preconditioned  iterative  method  to  approximate  from  (3.2). 

As  an  initial  guess  for  for  n  ^  2,  we  shall  use  quadratic 

extrapolation.  Specifically,  we  shall  use 


(3.8) 


«  n  -  n-1  .  n-2 

3£q  *  2y  -  3y  +  Y 


as  the  initialization  for  the  iterative  procedure  for  Y*'^^  “  Y*'*  Since  we 

n  n-1  .  n-2  .  .  ^  ^  n+1 

use  Y  »  Y  »  and  Y  rn  the  coefficient  matrices  to  determine  y  , 

the  errors  in  the  approximate  solution  will  accumulate. 

In  order  to  estimate  the  cumulative  error,  we  first  consider  the 


single  step  error.  Define  Y  to  satisfy 


(3.9)  l'‘(y)(y"'^^ 


-  y")  -  [< 


c  .  iMli  A-,,,] 


-  Y  )  *  I^(Y)  *  n  1  , 


froo  (3.2) .  For  all  of  the  analysis  to  follow,  we  can  use  any  precondi¬ 
tioned  iterative  method  which  yields  norm  reductions  of  the  form 


(3.10,  i0J|o^Y>^Y"‘  -  3y"  *  3y"‘’  -  Y-'^  II,  . 

where  0  <  <  1  and  the  subscript  iniJicates  the  Euclidean  norm  of  the 

vector.  A  particularly  efficient  iterative  procedure  for  obtaining  (3.10) 
is  the  preconditioned  conjugate  gradient  method  presented  in  11,2,9,10,13]. 


(3.11) 


a)  Ikll^  -  (cV'V>) 

b)  a(v")Vv>,7v5) 


c)  Ill^lIL  =  Ikll^  n  ' 

"  a 

be  special  norms  and  semi-norms.  Note  that  |1  •  H  ^  are  uniformly 

a 

equivalent  to  ll^*l|*  Then  letting 

M 

(3.12)  ^  -  I  » 

i-1 


(3.12) 

ve  see  that 

yn+l 

f  ^ 

c  — 

(3.13) 

1 

(f(t",o”),x)  +  <g(t“),x>. 


fc^n.^  0-1  1  ^  r 

(3.13)  ^  (At)  J  (  i 

-  (f(t",o”),x)  +  <g(t“),x>.  X  «  * 

Ke  also  see  that,  using  (3.11),  (3.10)  can  be  written  as 

(3.14)  1I|V”*'  -  V”**|l|„  i  Pilll«V|ll„.  n  .  2  , 


where 


(3.15) 


*> 

...  ,  n  _  n+1  n 

b)  s  ^  ~  <fi  t 

c) 

,  5  -  3,”  .  3,"-*  -  V."-’ 
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0  1  2 

Ke  next  discuss  a  starting  procedure  for  obtaining  V  ,  V  ,  and  V  . 

0  J 

Ke  shall  follow  the  ideas  of  Ill]  in  determining  V  and  V  .  Let 

=  W(0);  i.e.  project  into  ll^.  This  will  require  the  factoriza¬ 

tion  of  one  additional  matrix  to  solve  tlie  elliptic  problem  (2.7)  .  Then 
approximate  u(x,At)  by 

u  =  u(x,0)  +  At  (x,0)  +  (x,0)  . 

^  3t‘' 

*  g2 

Project  u  into  M.  .  The  derivative  — ^  is  evaluated  using  the 

^  at^ 

differential  equation.  The  solution  of  the  second  elliptic  problem  can  be 
approximated  using  the  factored  matrix  used  to  determine  V^.  We  can 


thus  obtain  the  estimate 


(3.16)  11V(V  -  W)°|l  +  llV(V  -  W)^ll  +  lld^(V  -  W)°||  <c{h^+  (At)^}  . 

0  1  2 

Once  we  have  V  and  V  satisfying  (3.16),  V  can  be  determined 
using  the  same  preconditioned  iterative  procedure  as  above  by  initializing 
the  iterative  procedure  by  *q  *  For  details  of  a  starting 

procedure  using  the  iterative  procedure,  see  (10]  . 

4.  A  Priori  Error  Estimates 

In  this  section  we  develop  a  priori  bounds  for  the  errors  u"  - 

and  V®  -  u"  for  the  procedures  defined  in  (2.12)  and  (3.13)  respectively. 

Similar  result.*  yielding  optimal  order  H^-estimates  can  be  obtained  using 

similar  techniques  for  the  procedures  defined  in  (2.13)  and  (2.14)  and  their 

2 

iterative  counterparts.  Theorem  4.1  yields  optimal  order  L  -estimates 
for  the  procedure  satisfying  (2.12)  and  (3.16)  under  restrictions  given 
in  (4.18).  Under  the  slightly  stronger  mesh-ratio  restriction 
(4.1)  At  <  C*h  , 


11- 


we  obtain  optimal  order  L^-cstimates  for  the  iterative  procedure  satisfy¬ 
ing  (3.13)  and  (3.16)  in  Theorem  4.2. 


Theorem  4.1.  Let  S,  Q,  R,  and  the  restrictions  on  of 

Section  2  hold.  Let  u”  satisfy  (2.12)  and  (3.16).  Then  there  exist 

constants  t»  h.,  and  K,  =  K-(K.;  i  =  0,...,5)  such  that  if  r  > 

wool  i. 

d/4 

At  ^  1,  h  £  h^,  and  At  <  h  , 

(4.2)  sup(||u  -  u|(  +  h||u-  1  +  (At)^}  . 


Proof .  Let  n”  “  u”  -  W**  and  C**  =  u”  -  W*'.  From  (1.1),  (2.7)  and 
(2.12),  we  see  that 

.n+1  _n-l  1  r 

L3t- 


(4.3) 


(cd^c".x)  +  [a{u")v  ^ '^xj  =  (°[^  ■  ‘*t”"]'x]  +  <n".x) 
+  |a(u*‘)Vw"  -  a(u")7  — - - -  ,7xj 


+  (f(t",u”)  -  f(t",u"),x)» 


X« 


We  shall  let  x  ■  ~  c”  ^  “  At(d^C**  +  ^)  in  (4.3).  Using  this 


test  function  and  (3.11),  the  left  hax\d  side  of  (4.3)  becomes 

,n+l  .  .n-1 


d_&  m  »n^l 

,.5  “  d.  c  \  f  _n+l  ^ 

[c  -  ,At(d^c"  +  d^C  )J  +  |a(o")7  ^ - 1- 


(4.4) 


-  IIv“Ilc  -  llv"'"llc  *1 


n-l||2  ^  1  rii_n+lif2 


a”  a 

Zn  order  to  obtain  telescoping  sums  when  (4.4)  is  summed  on  n,  we  must 
shift  the  indices  in  two  of  the  terms  above.  Mote  that 

.n-lii2  ^  ii_a-l||2  ^  __n-l. 


llv  11  n  *  II'  H  n-2  ”  *W**”*)m 


(4.5)  -  ^lu 
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Wc  shall  henceforth  use  C  as  a  generic  constant  in  our  analysis.  For 
the  first  term  on  the  right  of  (4.3) »  we  obtain 

if-1  !  _  _ ,  ■x  I 


(4.6) 


We  next  bound  the  second  and  fourth  terms  on  the  right  of  (4.3)  as  follows 


I  (n”  +  f(t",u")  -  f (t”,u") ,At(d^;"  +  d^?"'^)) 
n=l 


(4.7) 


n»*0 

We  split  the  third  term  on  the  right  side  of  (4.3)  as  follows 


^  w"-^l 


(4.8) 


I  L(u")7fw"  -  5?— ^ 
1-1 

•  I  i  tTj  +  1  • 


n  n  w""^^  + 

+  [a(u”)  -  a(U  )]7  - - - 


.vx)l 


In  order  to  treat  the  terms  in  (4.8)  t  we  shall  sum  by  parts  in  time, 

I T  ,7  C  +  c")  -  (c"  +  c"“^)  J )  I 

•n-l  ' 

<  |Y  -  a(u"“^)17|w"  -  ].7(c"  +  c“'^)] 

.  |Y  -  {^>-£4^1]..,.' 


(4.9) 


c-,]l 

.}].„(%  c-.]l 


♦  |(a(u‘“^7(|  «V”S,7(C^-tC^“S)|  +  |(a(u^)7(|  «V).7(^^  +  C°))| 


-®{ji  *  “'’1 


♦itIlc'llVi  *  Ilc‘-‘IIV2  ♦  ll'‘ll%  *  • 


.t-ln2 


13 


similarly,  we  sec  that 
t-1 


(4.10) 


'n=2 


|i  -  (c” . 

•n=l  ' 

|la(u")  -  a(u")  -  {3(0"“^  -  a(u"'S}lV^^ - ^ - .  V(i;"  +  c"'^  ]  | 

4  j|(a(u^"^)  -  a(U^"^)17  - ,  ^)j| 

♦  j|(a(u^)  -  a(U^)]V  ”  2  1 

Combining  (4.3)- (4.10) ,  we  see  that  after  summing  (4.3)  on  n  for  n  ®  1 
to  n  «  1  -  If  we  use  Lemma  2.1  to  obtain 

A  A 

n^O  L  '  a  a 

(4.11) 

*  c  Y  (iiv"iic  *  ii'"iiYi  *  ii'”iic}“ 

n-O  '  a  ^ 

♦  c,{||c‘||’  .  ||«‘-'||^}  *  c{||t'’|l|o  "  ‘“‘1 

Zn  order  to  bound  the  terms  multiplied  by  Cj^  on  the  right  side  of  (4.11) 
and  to  introduce  an  L^  term  on  the  left  hand  side  of  (4.11)  #  we  note  tl.at 

IU"*'ll’  -  llc"ll’  -  2At(ca  c-.c")  ♦  (At)’||d^;"|l’ 

(4.12) 


14 


Sum  this  iti<.quality  from  n  =  1  to  the  upper  limits  £  -  1  and  Z  -  2; 
then  multiply  the  resulting  inequalities  by  add  them  to  (4.11) 
and  use  (3.16)  to  obtain 


(4.13) 


'a  a  ' 

n=l  L  **  a  a  ■' 


2r  4 

h  +  (At) 


In  order  to  apply  the  discrete  Gronv/all  lemma  to  (4.13),  we  wish  to  show 

that  there  exists  a  constant  >  0  such  that 

0 


(4.14) 


T  .  1  . 

n*0  L 


The  given  starting  procedure  yields 


(4.15) 


.IS  • 

L 


He  shall  use  an  induction  argument  as  in  [24,10,13]  to  yield  (4.14)  with  the 
summation  starting  at  n  =  1.  For  1=2,  the  inequality  (4.13)  and  the 
estimate  (3.16)  imply  that 


(4.16)  *  * 

<  C{h^'  +  (At)^}  . 

Then  we  have  by  (2. 2. a),  (4.15),  and  (4.16), 


(4.17) 


ll«c^ll  .  *  llsc®!!  .  *  S 

L  I. 

<  C&th  ^ (h'  +  (At)^}  +  Cj  . 
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Then  if 


Then  if  (4.18)  is  satisfied  and  At  and  h  are  sufficiently  small «  our 
induction  argument  is  completed.  Then  since  (4.21)  holds  for  1  ^  1  ^  Nr 
using  (4.21)  r  (2.8)  and  the  triangle  inequality,  we  obtain  the  desired 
result  (4.2).  // 

We  shall  next  obtain  the  same  order  asymptotic  error  estimates  as 
derived  in  Theorem  4.1  for  the  approximation  V  defined  in  Section  3.  We 
shall  see  in  Section  5  that  the  worK  estimates  for  the  approximation  V  are 
far  superior  to  those  for  the  approximation  U  analyzed  above. 
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Thcorctn  4.2.  Let  S,  Q,  R,  and  the  restrictions  on  (Af^}  of 
Section  2  hold.  Let  v"  satisfy  (3.13),  (3.1G)  and  (3.14)  where 


■i  i  *  -ilr] }  “ 


V 

Then  there  exist  constants  t,  h^^,  and  (C  ,K^,  i  *  0,...,5)  such 

that  if  r  >  “,  At  ^  T,  h  £  h^,  and  At  ^Tnin{h^^^,C  h}, 

(4.22)  sup{(|u  -  V||  +  h||u-  1  K^{h’^  +  (At)^}  . 

t" 

Proof:  Let  z”  =  v”  -  w”  and  n"  be  as  eibove.  From  (1.1),  (2.7), 


and  (3.13),  we  see  that 


(cdV,x)  +  |a(v")7^ - - ,VxJ 

+  |a  (u")  Vw"  -  a  (V")  7  2  ^  ^  +  (f(t”,v")  -  f(t",u"),x) 

+  fc  »x]  +  (a(v")  |7(z"'^^  -  z"‘^^),7x),  X  «  • 


(4.23) 


We  note  that  except  for  the  last  two  terms  on  the  right  of  (4.23),  equation 
(4.23)  corresponds  exactly  with  (4.3).  We  must  thus  only  show  how  the 
last  two  terms  on  the  right  of  (4.23)  are  bounded.  From  (4.1),  (3.11) 


and  (3.14)  we  see  that 


4 


.  i  III,-:  - 


(4.24) 


<Pl  [’ 


•  * 

a  KqC 


2c 


*  2 

*  V 


^][l|a,i"ll,  *  lU,t''-MlJllUV|||„ 


1^]  [ll^t'llc*  IIV''"llc][llV''llc*  '1^'"“". 


Il'lt^"'’llc  ♦  ClitI*]  • 


We  then  see  that  if 


(4.25) 


•  *2-1 
a  K  C 


then 

(4*26)  I  I  ||d^;”||^  +  C(At)^  . 

n»l  n“0 

The  rest  of  the  proof  follows  as  in  the  proof  of  Theorem  4.1.  // 

Similar  techniques  can  be  used  to  give  a-priori  error  estimates  for 
the  approximations  given  in  (2.13)  and  (2.14)  as  well  as  for  the  correspond¬ 
ing  iterative  approximations  defined  in  Section  3.  Since  in  the  major 
applications,  the  coefficients  depend  upon  Vu  (the  strain),  the  techniques 

presented  here  will  only  yield  optimal  order  H^-estiaates  instead  of  the 
2 

optimal  order  L  -estimates  obtained  in  Theorem  4.1  and  4.2. 


5.  Computational  Considerations 

Zn  this  section  we  shall  consider  some  rough  operation  counts  to 
estimate  the  computational  complexity  of  the  methods  presented  here.  We 
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shall  show  that  the  iterative  methods  presented  in  Section  3  allow  us  to 
obtain  near  optimal  order  work  estimates.  Therefore,  those  methods  are 
very  efficient  computationally. 

First  consider  d  =  2  and  the  second  order  hyperbolic  equation.  Let 
M  be  the  dimension  of  and  N  be  the  number  of  time  steps.  George  [17] 

has  shown  in  some  special  cases  that  the  procedure  of  setting  up  and 
factoring  l"  (from  (3.9))  requires  0(M^^^)  operations  and  that  the  solu¬ 
tion  of  (3.2),  given  the  factorization,  requires  0(M  log  M)  operations. 
Hoffman,  Martin  and  Rose  [19]  have  shown  that  such  bounds  are  minimal. 
n>erefore,  if  we  conjecture  the  validity  of  the  above  estimates  for  our 
problem  and  refactor  l"  and  solve  (3.2)  at  each  time  step,  the  total  amount 
of  work  done  is 

(5.1)  +  M  log  M})  =  0(NM^'^^)  . 

We  note  that  the  work  of  factorization  dominates  the  estimates. 

Using  the  preconditioned  iterative  methods  presented  in  Section  3,  one 

does  not  have  to  refactor  at  every  time  step.  With  d  ■  2  we  have 

_  r  £ 

(5.2)  M  (At)”^  »  h  ^  . 

We  are  willing  to  refactor  periodically,  but  our  goal  is  to  have  the  total 

work  estimate  (5.1)  dominated  by  the  work  of  solving,  0(NM  log  M) .  We  shall 

•««  that  for  r  ^  3,  this  goal  can  be  achieved.  For  r  *  2  (piecewise 

linear  elements),  the  work  of  factoring  one  matrix  is  already  almost  as 

3/2 

large  as  the  total  work  of  solving  (0(M  log  M)  in  this  case).  If  we 
refactor  and  update  the  preconditioning  matrix  N  equally  spaced  times, 
one  can  show  that  the  preconditioning  matrices  are  sufficiently  comparable 
to  the  true  matrices  that  each  iteration  of  the  iterative  procedure  yields 


a  norm  reduction  of  0(  (At) .  Then  for  At  sufficiently  s-Hvill,  five 
iterations  per  tine  stop  will  satisfy  the  norm  reduction  reciuirenont  of 
(4.25).  Next,  for  n  *  3  (piecewise  quadratic  elements),  (5.1)  and  (5.2) 
show  that  the  total  wor)c  is 

1  3  1  1 

(5.3)  0(N^M^  +  5NM  log  M)  =  0(M^^  ^  +  m'*  log  M)  =  0(M^  log  M)  , 

and  the  work  of  solving  dominates  the  estimate.  If  r  =  2,  the  total  work  is 

1+1  3  12 

(5.4)  0(M^  ®  +  5M^  log  M)  =  0(M®)  , 

n 

which  IS  still  much  better  than  the  0(M  )  work  estimate  if  the  matrices 
are  factored  at  each  time  step. 

If  i  4,  one  can  refactor  and  update  the  preconditioning  matrix  more 

1/2 

frequently  (specifically  N  equally  spaced  times),  obtain  a  norm  reduc- 
1/2 

tion  of  0((At)  )  with  each  iteration  and  by  iterating  only  three  times 

per  time  step,  still  have  the  work  of  solving  dominate  the  work  estimate. 

If  r  «  4  (piecewise  cubics) ,  the  total  work  is 

7 

(5.5)  0(M^  +  3M^  log  M)  «  0(M^  log  M)  . 

We  thus  see  that  If  r  2  3,  then  by  refactoring  and  updating  the  pre¬ 
conditioning  matrix  sufficiently  often  (depending  upon  r)  ,  the  total  work 
is  of  the  order  0(NM  log  M) .  Since  the  total  number  of  unknowns  in  the 
problem  is  0(NM),  we  see  that  we  can  obtain  almost  optimal  order  work 
estimates  for  r  2 

2 

For  d  *  3,  the  work  of  factoring  a  matrix  is  0(M  )  while  the  work  of 
solving  the  result  is  0(M^^^).  Thus  if  1^2  the  total  work  of  solving 
again  dominates  the  work  of  factoring  a  matrix.  Thus  if  refactoring  is  done 
sufficiently  infrequently  (depending  upon  r)  the  total  work  of  solving 
will  again  dominate  the  total  work  estimates. 
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It  is  computationally  wasteful  to  iterate  sufficiently  na:iy  timi.-c  at 
each  time  seep  to  achieve  the  pessimistic  bounds  given  by  (4.25).  Instead, 


one  can  monitor  the  norm  reduction  actually  produced  at  each  step  of  the 
Iteration  and  stop  iterating  when  sufficient  norm  reduction  is  achieved. 
Additional  stopping  criteria  can  be  imposed  in  this  monotoring  process. 
[10)  for  a  discussion  of  stepping  criteria  for  a  related  problem. 


See 


-21- 


KEi'ERi'-NCL'S 


1.  O.  Axelsson,  "On  preccr.ditioning  and  convergence  acceleration  in  sparse 
matrix  prcblens,"  CER.N'  European  Organization  for  Nuclear  Research, 

Geneva,  1974. 

2.  O.  Axelsson,  "On  the  computational  complexity  of  some  matrix  iterative 
algorithms,"  Report  74.06,  Dept,  of  Computer  Science,  Chalmers  University 
of  Technology,  Goteberg,  1974. 

3.  G.  A.  Baker,  "Error  estimates  for  finite  element  methods  for  second 
order  hyperbolic  equations,"  SIAM  J.  Numer.  Anal,  13  (1976),  pp.  564-576. 

4.  J.  H.  Bramble  and  P.  H.  Sammon,  "Efficient  higher  order  single-step 
methods  for  parabolic  problems,"  (to  appear). 

5.  B.  D.  Coleman  and  W.  Noll,  "The  thermodynamics  of  elastic  materials  with 
heat  conduction  and  viscosity,"  Arch.  Rat.  Mech.  Anal.  13  (1963), 

pp.  167-178. 

6.  C.  M.  Dafermos,  "The  mixed  initial-boundary  problem  for  equations  of  non¬ 
linear  viscoelasticity,"  J.  of  Diff.  Eqn.  6  (1969),  pp.  71-81. 

7.  J.  C.  Dendy,  Jr.,  "An  analysis  of  some  Galerkin  schemes  for  the  solution 
of  nonlinear  time-dependent  problems,"  SIAM  J.  Kumer.  Anal.  12  (1975), 
pp.  541-565. 

8.  J.  E.  Dendy,  Jr.,  and  G.  Fairweathcr,  "Alternating-direction  Galerkin 
methods  for  parabolic  and  hyperbolic  problems  on  rectangular  ivolygons," 
SIAM  J.  Numer.  Anal.  12  (1975),  pp.  144-163. 

9.  J.  Douglas,  Jr.,  and  T.  Dupont,  "Preconditioned  conjugate  gradient 
iteration  applied  to  Galerkin  methods  for  a  mildly  nonlinear  Dirichlct 
problem,"  Sparse  Matrix  Calculations ,  J.  R.  Bunch  and  D.  J.  Rose,  cds.. 
Academic  Press,  New  York,  1976,  pp.  333-348. 


-22- 


rassHw-c 


10.  J.  Douglas,  Jr.,  T.  Dupont,  and  R.  E.  E'-^ing,  "Incomplete  iteration  for 
time-stepping  a  Galcrkin  method  for  a  quasilincar  parabolic  problem," 

Sl^o;  J.  N'u.mer.  Anal,  (to  appear). 

11.  T.  Dupont,  "L^-estimates  for  Galerkin  methods  for  second  order  hyperbolic 
equations,"  SI.^!  J.  Numer.  Anal.  10  (1973),  pp.  880-889. 

12.  T.  D-jpont,  G.  Fairveather,  and  J.  P.  Johnson,  "Three-level  Galerkin 
methods  for  parabolic  equations,"  SIAM  J.  Numer.  Anal.  11  (1974), 
pp.  392-410. 

13.  R.  E.  Ewing,  "Time- stepping  Galerkin  methods  for  nonlinear  Sobolev  partial 
differential  eq’uations,"  SIAM  J.  Numer.  Anal.  15  (1978),  pp.  1125-1150. 

14.  R.  E.  Ewing,  "Efficient  time-stepping  procedures  for  miscible  displacement 
problems  in  porous  media,"  Mathematics  Research  Center  Technical  Summary 
Report  #1934,  University  of  Wisconsin-Madison,  Madison,  Wisconsin. 

15.  R.  E.  Ewing,  "Efficient  time-stepping  methods  for  miscible  displacement 
problems  with  nonlinear  boundary  conditions,"  (to  appear). 

16.  R.  E.  Ewing  and  M.  F.  Wheeler,  "Galerkin  methods  for  miscible  displacement 
problems  in  porous  media,"  Mathematics  Research  Center  Technical  Summary 
Report  #1932,  University  of  Wisconsin-Madison,  Madison,  Wisconsin. 

17.  A.  George,  "Nested  dissection  on  a  regular  finite  element  mesh,"  SIAM  J. 
Nhimer.  Anal.  10  (1973),  pp.  345-363. 

18.  J.  M.  Greenberg,  U.  C.  MacCamy,  and  V.  J.  Mizcl ,  "On  the  existence, 

uniqueness,  and  stability  of  solutions  of  tlio  equation  o' (u^)u^^^  +  ~ 

p^u^^,"  J.  Math,  and  Mech.  17  (1963),  pp.  707-728. 

19.  A.  J.  Hoffman,  M.  S.  Martin,  and  D.  J.  Rose,  "Complexity  bounds  for 
regular  finite  difference  and  finite  element  grids,"  SIAM  J.  Numer.  Anal. 

10  (1973) ,  pp.  364-369. 


-23- 


20.  M.  Lighthill,  "Dynaraics  of  rotating  fluids:  a  survey,"  J.  Fluid  Mech. 

26  (1966),  pp.  411-436. 

21.  A.  E.  H.  Love,  A  Treatise  on  the  Mathematical  Theory  of  Elastic! ty , 

Dover,  New  York,  1944. 

22.  R.  C.  MacCamy,  "Existence,  uniqueness,  and  stability  of  solutions  of  the 

3 

equation  u  =  —  (o(u  )  +  1 (u  )u  .),"  Report  68-18,  Dept,  of  Math., 

ww  oX  X  X  Xu 

Carnegie  Inst,  of  Technology,  Carnegie-Mellon  University. 

23.  G.  W.  Platzman,  "The  eigenvalues  of  Laplace's  tidal  equations,"  Quart. 

J.  of  Royal  Meteorological  Soc.  94  (1968) ,  pp.  225-248. 

24.  H.  H.  Rachford,  Jr.,  "Two-level  discrete-time  Galerkin  approximations  for 
second  order  nonlinear  parabolic  partial  differential  equations,"  SIAM 

J.  Numer.  Anal.  10  (1973),  pp.  1010-1026. 

25.  R.  E.  Showalter,  "Regularization  and  approximation  of  second  order 

evolution  equations,"  SIAM  J.  Math.  Anal.  7  (1976),  pp.  461-472. 

2 

26.  M.  F.  Vfheeler,  "A  priori  L  -error  estimates  for  Galerkin  approximations 
to  parabolic  partial  differential  equations,"  SIAM  J.  Numer.  Anal.  10 
(1973),  pp.  723-759. 


REE/ j vs 


-24- 


security  classification  of  this  pace  (Whmn  Dmim  Bnl»r»dJ 


REPORT  DOCUMENTATION  PAGE 


r  report  number  12  govt  accession  no. 

#1996  _  , _ 


4.  TITLE  (tmd  Sublllld) 

ON  EFFICIENT  TIME-STEPPING  METHODS  FOR  NONLINEAR 
SECOND  ORDER  HYPERBOLIC  PARTIAL  DIFFERENTIAL 
EQUATIONS 


7.  AUTMORr«J 


Richard  E.  Ewing 


».  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Mathematics  Research  Center,  University  of 

610  Walnut  Street  Wisconsin 

Madison.  Wisconsin  53706 _ 


11.  controlling  office  name  ano  address 


See  Item  18. 


I 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


3  RECIPIENT’S  CATALOG  NUMBER 


5.  TYPE  OF  REPORT  A  PERIOD  COVERED 

no  specific 
reporting  period 


6  PERFORMING  ORO.  REPORT  NUMBER 


8.  CONTRACT  OR  GRANT  NUMBERCl/ 

MCS78-09525 
DAAG29-7  5-C-0024 
DAAG29-78-G-0161 


10.  PROGRAM  element.  PROJECT.  TASK 
AREA  A  WORK  UNIT  NUMBERS 

Work  Unit  Number  7  - 
Numerical  Analysis 


12.  report  date 
September  1979 


13.  NUMBER  OF  PAGES 

24 


MONITORING  VGCNCY  NAMC  A  dtUmrmit  trom  Controlling  Ollier)  'S.  SECURITY  CLASS,  (of  thio  roport) 

UNCLASSIFIED 


IS«.  DECLASSIFICATION  DOWNGRADING 
SCHEDULE 


16.  DISTR  euTlON  STATEMENT  (ol  thia  Roport) 

Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (el  the  abmtrect  entered  In  Block  20,  II  dllleteni  from  Report) 


1$.  SUPPLEMENTARY  NOTES 

U.  S.  Army  Research  Office  and  Nationa; 

P.  O.  Box  12211  Washing! 

Research  Triangle  Park 
North  Carolina  27709 


19.  KEY  WORDS  (Continuo  on  rovmrao  aido  if  nacaaaary  and  idantify  by  block  numbar) 

Galerkin  methods,  error  estimates,  hyperbolic  equations. 


National  Sci 'nee  Foundation 
Washington,  D.  C.  20550 


20.*^OSTRACT  (Contfnua  on  rararaa  aida  II  nacaaaary  and  Idantify  by  block  numbar) 

techniques  useful  for  efficiently  time-stepping  Galerkin  methods  for  various 
types  of  time-dependent  partial  differential  equations  are  presented  and 
analyzed.  Second-order  quasilinear  hyperbolic  problems  with  smooth  solutions  are 
studied  as  a  simple  model  problem  for  illustrating  the  widely  applicable  tecli- 
niques.  The  procedure  involves  the  use  of  a  preconditioned  iterative  method  for 
approximately  solving  the  different  linear  systems  of  equations  arising  at  each 
time-step  in  a  discrete-time  Galerkin  method.  Optimal  order  spatial  errors 
and  almost  optimal  order  work  estimates  are  obtained  for  the  second-order  hyper- 


FORM 
I  JAN  73 


EDITION  OF  1  NOV  •>  IS  OBSOLETE 


...  _  bolic  equation 

UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Dele  Entered) 


